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Abstrac  *. 


Let  G=(V,E)  denote  an  undirected  network  with  node  set  V  and 
arc  set  E={1,...,N}  .  Arcs  fail  randomly  and  independently  with 
probability  1-q^  for  V  i«E  .  This  paper  describes  a  numerical 
procedure  for  estimating  g(s,t)  ,  the  probability  that  nodes  s  and 
t  are  connected  for  specified  s,teE  ,  with  bounds  on  absolute  error 
proportional  to  1/K  for  a  specified  nonrandom  finite  sequence  and 
proportional  to  log  K/K  for  certain  nonrandom  infinite  sequences, 
where  K  is  the  number  of  replications.  These  convergence  rates  are 
best  possible.  Although  the  infinite  sequences  have  a  slower  convergence 
rate  than  the  finite  sequence  has,  they  offer  the  convenience  of 
allowing  one  to  add  replications  and  retain  the  rate  whereas  the  finite 
sequence  does  not.  These  bounds  improve  on  the  convergence  rate 
OO/K1^2)  for  the  standard  error  in  the  case  of  independent  Monte  Carlo 
replications  based  on  random  sampling.  Moreover,  they  hold  with 
certainty.  Algorithms  for  computing  estimates  are  shown  to  have  an 
upper  bound  0(N)+0(max(N, | V| ))  on  time  complexity  per  replication 
as  K  +  »  . 

The  paper  first  describes  the  estimation  of  g(s,t)  for  q1=...=qN= 
by  using  a  tabled  truncated  binomial  distribution  together  with  the 
A-canonical  representation  of  a  positive  integer.  It  then  describes 
how  to  incorporate  stratified  sampling  to  estimate  g(s,t)  as  a 
function  of  p  at  small  marginal  increase  in  time  complexity.  Next, 
the  paper  extends  the  stratified  sampling  method  to  the  case  of 
unequal  qlt...,qN  . 

An  example  based  on  a  network  of  30  arcs  illustrates  the  techniques. 
Lastly,  the  paper  extends  the  method  to  the  estimation  of  g(s,T)  ,  the 
probability  that  s  and  t  are  connected  for  all  tcTCE-(s)  . 

KEY  WORDS:  Monte  Carlo  methods,  network  connectedness, 
network  reliability,  quasi  random  points. 


Introduction 


Let  G=(V,E)  denote  an  undirected  network  with  node  set  V 
and  arc  set  E={1,...,N)  .  Arcs  are  assumed  to  fail  randomly  and 
independently  with  probability  1-q^  for  all  ieE  .  We  wish  to 
estimate 

g(s,t)  =  the  probability  that  nodes  s  and  t 
are  connected  for  specified  s,teV  . 

It  is  well  known  that  direct  evaluation  of  g(s,t)  takes  0(2 IEl) 

steps,  a  computation  that  can  be  infeasible  even  with  a  moderate 

number  of  arcs.  At  least  three  approaches  have  been  suggested  to 

reduce  the  severity  of  this  difficulty.  The  first  concentrates  on 

networks  with  special  structure.  For  example,  see  Rosenthal  (1977). 

The  second  relies  on  finding  bounding  inequalities.  For  example,  see 

Ball  and  Provan  (1983)  and  Zemel  (1982). 

The  third  uses  estimation  techniques  among  which  the  most  common 
are  Monte  Carlo  procedures  based  on  independent  trials  or  replications 
that  use  random  sampling.  Van  Slyke  and  Frank  (1972)  describe  how 
Monte  Carlo  methods  using  a  stratified  sampling  design  lead  to  estimators 
of  network  reliability  with  considerably  smaller  variances  than  would 
obtain  if  pure  random  sampling  were  employed.  Easton  and  Wong  (1980) 
demonstrate  a  similar  benefit  using  conditional  sampling  and  Kumamoto, 
Tanaka  and  Inoue  (1977)  show  a  like  benefit  for  dagger  sampling,  a  special 
form  of  antithetic  variates  (Hammersley  and  Handscomb  1964).  Karp  and 
Luby  (1983)  describe  a  Monte  Carlo  procedure  based  on  importance 
sampling  whose  computation  time  complexity  to  achieve  a  specified 
accuracy  is  linear  in  the  number  of  cutsets. 


»’ArA7» 


While  each  exploits  network  structure  in  a  specialized  way,  these 

Monte  Carlo  procedures  share  one  shortcoming  in  common.  They  all  lead 

1  /2 

to  standard  errors  that  converge  as  0(1/K  '  )  where  K  is  the  number 
of  independent  trials  or  replications.  To  improve  this  convergence 


rate  Fishman  (1983)  describes  a  numerical  procedure  that  uses  quasi  random 
point  sequences  together  with  a  conditionality  property  based  on  cutsets 
to  produce  an  estimate  of  g(s,t)  whose  absolute  error  has  an  upper 
bound  proportional  to  (log  K) 1  '/K  where  P  denotes  the  union  of 

all  minimal  s-t  pathsets  and  R  is  a  minimal  cutset  with  the  property 
that  |R ns |=1  for  each  s-t  pathset  S  .  Moreover,  since  these 
sequences  are  nonrandom,  the  bound  holds  with  certainty.  These  results 
imply  that  for  a  given  network  6  quasi  random  point  sequences  produce 
a  deterministic  upper  bound  on  absolute  error  that  converges  faster 


than  1/K 


1/2 


The  purpose  of  this  paper  is  to  describe  related  estimation  techniques 


that  improve  4h+*  convergence  rate  to  1/K  when  a  specified  nonrandom 
finite  sequence  of  sample  points  is  used  and  to  (log  K)/K  when 
certain  nonrandom  infinite  sequences  of  sample  points  are  used.  In 
addition,  algorithms  are  given  for  computing  estimates  of  g(s,t)  with 
1/K  and  (log  K)/K  convergence  rates  with  computation  time  complexities 
per  replication  having  an  upper  bound  0(N)  +  0(max(N, I V 1 ) )  as  K  -*■  <*>  . 


Although  the  infinite  sequences  have  a  slower  convergence  rate  than 


the  finite  sequence  has,  they  offer  the  convenience  of  allowing  one 
to  add  replications  as  desired  and  retaining  the  rate.  By  contrast, 
the  finite  sequence  does  not  allow  this  addition,  once  K  is  Initially 


Section  1  describes  the  technique  for  the  special  case  of 
equal  failure  probabilities.  In  particular,  it  shows  how  sampling 
from  a  tabled  truncated  binomial  distribution  and  applying  the 
A-canonical  representation  of  an  integer  B  (Kruskal  1963  and 
Katona  1966)  make  these  convergence  rates  and  time  complexity  possible. 
Section  2  then  shows  how  to  incorporate  stratified  sampling  to  allow 
the  estimation  of  g(s,t)  as  a  function  of  p  at  relatively  small 
increase  in  time  complexity  while  retaining  the  desirable  convergence 
rates.  Section  3  shows  how  the  method  of  stratified  sampling 
facilitates  the  estimation  of  g(s,t)  for  the  case  of  unequal 
q^,...,qN  ,  while  preserving  the  convergence  rates  and  time  complexity. 
Section  4  illustrates  the  technique  of  Section  1  for  a  network  of 
30  arcs.  Lastly,  Section  5  extends  the  analysis  to  the  case  of 
estimating  g(s,T)  ,  the  probability  that  s  is  connected  to  all 


V 
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1.  The  Basic  Method 


For  all  nodes  u,veV  define: 

C(u,v)  =  minimal  u-v  cutset  of  minimal  cardinality 
M(u,v)  =  minimal  u-v  pathset  of  minimal  cardinality 
L(u,v)  =  set  of  all  minimal  u-v  pathsets 
P(u,v)  =  union  of  all  minimal  u-v  pathsets 


and 


V(u,v)  =  node  set  corresponding  to  P(u,v)  . 


Let 


Fjfn.e)  •  1  (")  eJ(l-e)"'j 

J  i=0  J 

n  =  integer  >  0  and  0  <  e  < 
0  s  j  s  n  . 

For  each  arc  1*1,... tN  let 

q.j  -  probability  that  arc  i  operates 


and 


2^  *  1  if  arc  i  operates 

*  0  if  arc  i  fails. 

Then  for  the  network,  let 

z  *  (z^,...,zN)  =  network  state 
Z  *  set  of  all  network  states  z 


1  -  set  of  all  network  states  with  z,+...+z„=m 
m  IN 


k(u,v|z)  *  1  -  H 


SeL(u.v) 


(1  -  n  Z.) 


leS 


Rju.v)  *  l  k(u,v|z) 

ZcZ 

m 


and 


R(u,v)  *  l_  k(u,v|z)  -  l  Rm(u,v)  • 


ZeZ 


m*l 


Note  that  given  z  k(u,vlz)  is  the  conditional  probability 
(0  or  1)  that  u  and  v  are  connected.  Also  Rm(u,v)  denotes 
the  number  of  distinct  states  z  that  induce  u-v  connectedness  when 
m  arcs  operate  and  R(u,v)  denotes  the  number  of  distinct  states  z 
that  induce  u-v  connectedness  regardless  of  the  number  of  operating 
arcs. 

Expression  (2)  is  of  seminal  importance  throughout  this  paper. 
Suppose  we  collect  data  x^,...,xK  on  K  sampling  experiments  or 
replications  where 

Xj  = 
and 

x. .  =  1  if  arc  i  operates  on  replication  j 

'  J 

*  0  otherwise. 

Here  x^,...,xK  are  sample  network  states  and 

gK(s,t)  =  X  ^  k(s ,t Ixj )  (5) 


is  the  proportion  of  trials  on  which  s  and  t  are  connected.  As  we 
show  shortly  the  sampling  plan  for  choosing  x^,...,xK  determines  how 
good  an  approximation  gK(s,t)  is  to  g(s,t)  . 

It  is  important  to  note  that  (2)  is  used  merely  for  representation 
and  analysis.  To  determine  the  value  of  k(u,vlx.)  on  each  replication, 

J 

we  assume  that  a  "depth-first  search"  algorithm,  as  in  Aho,  Hopcroft 
and  Ullman  (1974,  pp.  176-177),  is  employed.  This  algorithm  has  an 
upper  bound  on  computation  time  complexity  0(lP(u,v)l,  I V(u,v) I )  where 


IS  I  denotes  the  cardinality  of  the  set  S  . 

Sampling  Sequence*  and  Exaoa.  Bound* 

We  now  turn  to  the  selection  of  a  sampling  scheme  for  generating 
the  data  x^,...,xK  .  In  particular,  we  would  like  a  scheme  that  enables 
us  to  estimate  g(s,t)  by  $^(s,t)  and  related  estimators  in  a 
computationally  efficient  manner.  To  make  an  informed  selection,  we 
first  need  to  consider  several  basic  ideas  from  the  theory  of 
equidistributed  points.  Let 

I[a.6)(x)  =  1  1f  “SX<S 

*  0  otherwise 


and  let  v1 ,  v2,...  denote  a  sequence  of  points  in  [0,1)  .  Then  the 
sequence  v^,  is  said  to  be  uniformly  distributed  or  equi distributed 

on  [0,1)  if  (Schmidt  1977,  p.  1) 


1  K 

lifl* v  I  lr  R\ 
K—  K  j=l  Lo.PJ 


<»j>* 


for  all  o  and  b  such  that  0so<bs1  •  Conversely,  if  v.j,v2,... 
are  equidistributed  in  [0,1)  then  (6)  holds.  As  a  measure  of  error 


in  estimating  B-a  one  has  the  discrepancy  measures 

1  K 

dk  s  sup  l  !r«  si  (vi>  '  ^_o)l 

K  0sa<Bsl  K  j=l  Lo’8;  3 


m 


m 


*  ^  K 

°K  =  sup  lK  ^  !rO  o^ViJ  *  ° 
K  Os  a  si  K  j=l  L°‘“'  3 


with  the  inequality  (Kuipers  and  Niederreiter  1974,  p.  91) 


Dk  s  Dk  s  2Dk  . 


Note  that  DK  and  DK  are  worst  case  measures  of  the  errors  of 
estimation. 

For  the  finite  sequence  ;  j=l,...,K)  it  is  known 


(Niederreiter  1978,  p.  972)  that 


Dk=1/K,  (9) 

which  is  best  possible.  In  this  paper  we  devise  a  sampling  scheme 
that  enables  us  to  write  §K(s,t)  and  related  estimates  as  linear 
combinations  of  indicator  functions  with  argument  v.  and  thus  are 

J 

able  to  realize  the  convergence  rate  1/K  . 

Although  one  can  achieve  this  rate  in  practice  for  §K(s,t)  in  (5), 
situations  can  arise  in  which  after  K  replications,  one  finds  that 
a  larger  sample  size  is  necessary.  Unfortunately  the  finite  sequence 
;  j*l , . . .  ,K)  does  not  allow  for  adding  points 

VK+1*  vK+2’* ' * *VK+J  and  ac*”ev1n9  a  bound  proportional  to  1/(K+J) 
on  discrepancy.  However,  it  is  known  that  there  exist  infinite  sequences 
VpV2»...  (Kuipers  and  Niederreiter  1974,  p.  125)  such  that 


DKsc  (log  K)/K 


a*' ■.  aw  .I'.'*.- - v  *>  •*> 


where  c  depends  only  on  the  sequence.  One  such  sequence  due  to 
van  der  Corput  (1935)  considers  the  unique  dyadic  expansion  of  the 
integer  n 

I  , 

n  =  l  a. 2  a.ctO.l)  0<i<I 

i=0  1  1 

where  Ll=  log2  nj  .  Then  the  radical  inverse  function 

♦,(">  *  I  a,2"1'’ 

c  i=l  1 

can  be  used  to  compute  the  van  der  Corput  sequence  (v  .=  <j>2  (n+j )  ; 
j=l,2,...}  for  which  unpublished  results  of  Tijdeman  (Niederreiter  1974, 
p.  973)  show  that  for  all  K^l 

°K  5  1o92k+1)/k  (11a) 

and 

lim  sup  (KD*  -  --  log2  K)  *  |  +  y  log2  3  .  (lib) 

Bejian  and  Faure  (1977)  show  that  (11a)  also  holds  for  DK  and  (11b)  holds 
with  equality  for  DK  . 

Schmidt  (1977,  p.  28)  shows  that  for  infinite  sequences  the  best 
possible  lower  bound  is 

Dk  *  0( (log  K)/K) 


so  that  the  van  der  Corput  sequence  achieves  the  fastest  convergence 
rate  with  respect  to  K  .  Moreover,  no  other  infinite  sequence  is 
known  that  has  uniformly  smaller  discrepancy  than  this  sequence  has 


m m  ■  "A  m  .*  n  ■■  .■.  *  ^rr^rrrr.  w  t*  ■m.w.v.'jv.  *  •*, 


i.' 
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(Kuipers  and  Niederreiter  1974,  p.  127).  Note  that  the  van  der  Corput 
sequence  is  one  form  of  a  quasi  random  point  sequence. 

It  is  important  to  note  that  both  the  finite  and  infinite  sequences 
discussed  here  are  nonrandom  so  that  the  equality  (9)  and  the  bound  (10) 
hold  with  certainty.  By  contrast,  independent  replications  using 
random  sampling  lead  to  (Chung  1949) 


lim  sup 


/2K  D„ 


=  1 


w.p.  1  , 


K  /log  log  K 

indicating  a  considerably  slower  rate  of  convergence.  Therefore, 
provided  that  we  can  effectively  control  the  time  complexity  of  our 
proposed  sampling  plans,  the  foregoing  finite  and  van  der  Corput 
sequences  should  be  preferred  to  independent  random  sampling. 
Canonical  Rzpm&cntalion 

As  we  show  shortly  the  sampling  method  that  we  adopt  requires 


that  we  be  able  to  use  a  single  number  v.  on  replication  j  to  select 

J 


the  network  state  z  from  Z  or,  as  in  Sections  2  and  3,  from  Z 


m 


To  make  this  selection  possible  in  a  relatively  efficient  manner 
we  rely  on  an  ext^mely  useful  mathematical  relationship.  Let  B  be 
a  positive  integer.  Then  there  exists  an  A-canonical  representation  of  B 
(Kruskal  1963  and  Katona  1 966 )+ 


B,  B,  ,  B 

b  *  (/)  +  (  £]>  +...♦(/> 


(12) 


where 


BA  ”  >8^**1 


rI  am  grateful  to  nqy  colleague.  Professor  Scott  Provan,  for  making  me 


aware  of  this  relationship. 


\  <* 


•  »  *  »  *»*  »  *  •  *  •  •  . 
t - »-« -» -i -V - 


Most  importantly,  this  representation  is  unique. 

To  determine  B^,  for  B>0  one  can  use  the  relationship 

(Ball  and  Provan  1983,  p.  256) 

B.  =  max[b:  (J)  <  B] 

A  A  (13a) 


B.  =  max[b:  (!;)  <  B  -  T 
1  1  j=T+l 


Bi 

(jJ)] 


i=A-l ,A-2,. . . , £ 


where 


A  B 

z  =  max[i :  B  =  £  (^)i  i*l] 

j-i  J 


(13b) 


Suppose  that  we  know  that  exactly  A  of  N  arcs  operate  on 

N 

replication  j  so  that  T  x. .  =  A  .  To  determine  exactly  which  of 

i=l 

the  N  arcs  operate  one  can  proceed  as  follows:  Sample  B  from  the 
discrete  uniform  distribution  on  {0,1,... ,(A)-1 }  .  The  quantity  B 
is  the  "number"  of  the  combination  of  A  operating  and  N-A  failed 
arcs  on  replication  j  .  Then  if  B=0  one  takes 


xi  j  *  ^  i=l ,. . .  ,A 
*  0  otherwi se . 


(14a) 


If  B>0  and  *=1  one  takes 


xij  - ' 


1«B+1 
m 

otherwise  . 


m=l ,. . .  ,A 


(14b) 


Lastly,  if  B>0  and  &>1  one  takes 


■ 1 


1*1,...,£-1  ,  B4+1,...,Ba+1 


otherwise. 


(14c) 


Note  that  these  assignments  are  unique  functions  of  A  and  B  . 

Equal  PAobabllitieJ) 

We  first  consider  the  case  of  s-t  connectedness  with  equal 
operating  probabilities  q^p  ieP(s,t)  .  For  convenience  of 
exposition  we  take  C=C(s,t)  ,  L=L(s,t)  ,  M=M(s,t),  P=P(s,t)  , 
and  Rm=Rm(s,t)  ,  and  assume  P=E  so  that  |P|=N  .  Here  the  number 
of  operating  arcs  has  the  binomial  distribution  with  probability 
distribution  function  (1)  with  n=N  and  e=p  . 

Algorithm  C.l  describes  how  to  conduct  the  K  replications  that 
lead  to  the  computation  of  §K(s,t)  as  an  estimate  of  g(s,t)  .  Note 
that  the  number  of  operating  arcs  A  is  determined  (step  2c)  from  a 
truncated  distribution.  To  appreciate  why  this  is  so  observe  that  if 
the  number  of  operating  arcs  on  a  replication  exceeds  N-|C|  ,  then 
s  and  t  are  connected  with  probability  1.  Conversely,  if  fewer  than  |M| 
arcs  operate  then  with  probability  1  s  and  t  are  not  connected.  Therefore, 
it  is  to  our  advantage  to  restrict  sampling  (step  2c)  to  the  uncertain 
connectedness  outcomes  when  the  number  of  operating  arcs  exceeds 
1 M | — 1  but  is  less  than  N-|C|+1  .  Also  note  that  A  can  be  determined 
on  each  replication  with  time  complexity  0(1)  using  a  procedure 
described  in  Fishman  and  Moore  (1981). 

Theorem  1  gives  the  implications  of  algorithm  C.l. 

Theorem  1.  Suppose  algorithm  C.l  is  used  to  compute  §K(s,t)  .  Then: 

(i)  For  the  finite  sequence  (v^  =  ;  j*l,...,K} 

lgK(s,t)  -  g(s,t)|  s  2HQ/K  (15a) 


where 
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n  |  v  |  ii  i  w  i  * 

Q  S  min{  l  Rm  ’  l  [( 

m=|M|  m  m=]M I  r 


>  -  U>  • 


05b) 


(ii)  There  exist  infinite  sequences  {v.  ;  j=l,2,...} 

J 

such  that  the  bound  in  (15a)  holds  with  c(log  K)/K 
replacing  1/K  and  where  c  deoends  only  on  the 
selected  sequence. 

(iii)  One  can  compute  each  iteration  of  step  2  with  an  uoper 
bound  0(Np)  +  0(N)  +  0(max(N, | V I ) )  on  time  complexity 
per  replication  as  K  -*•  •. 

Proof.  For  (i)  observe  that  the  number  of  operating  arcs  A  has  a 

truncated  binomial  distribution  (step  2c)  on  the  integers  |M] ... . ,N- |C| 

with  parameters  N  and  p  .  Also,  for  given  A  it  is  clear  that  any 
N 

of  the  (A)  combinations  of  A  operating  arcs  are  equally  likely 
(step  2d).  The  quantity  B  merely  denotes  the  number  of  the  combination 
(given  A)  selected. 

Let  m*A  and  suppose  that  the  (JJJ)  possible  combinations  of 

operating  arcs  are  numbered  1 ,...,( Jj)  .  Let  JM  ■  (1 ,...,( jj)}  denote 

the  set  of  all  these  indices,  let  J*  denote  the  subset  of  J  for 

m  m 

which  s  and  t  are  connected  and1,  finally,  let  F  =F„(N,p)  .  Then 

mm 

the  unique  A-canonical  representation  of  B  leads  to 


N-ICI 


Ms.tlxJ  *  l  l  *  Irt-1  l  x  (  — 1 1 - 

m  'm 


Hvj”Fm-l+F|M| -1 


N-)CI 

*  I  1  *  ^ra  k  )  (vj) 
m*|M|  t€J  Latm  *  D*m'  J 


.-.v.  ^VV’AW.-V 


Algorithm  C.1 

Purpose:  To  compute  g  (s,t)  as  an  estimate  of  g(s,t)  . 

I\ 

Given:  Nt  p,  ICI  ,  IMI  ,  K  ,  s  and  t. 

1.  Set  parameters:  n=N,  e=p  and  S=0. 

2.  Perform  K  replications: 

a.  For  1=1,. ..,n  set  x..=0. 

*  J 

b.  Select  v.  . 

J 

c.  Determine  the  number  of  operating  arcs: 

A=min{m:  [Fm(n,e)-F|H|_1(n,e)]/H>v^.;  |M|<m<n-|C|} 

where  H  s  Fn.|c|(n*e^FlMl-l^n*6^  * 

d.  Determine  the  "number"  of  the  combination  of  A  operating 
and  N-A  failed  arcs: 


B  = 


L  9A(l-9)"-A  j 


e.  Determine  the  A-canonical  representation  t  ,  B£,...,BA  of 
using  (13).  If  B=0  ,  Jt=A+l  . 

f.  If  fc>l  ,  then  for  v=l  1  ;  By=v-1  . 

g.  Set  operating  indicators: 

For  v*1,...,A  ;  set  isBy+l  and  x^=l  . 

h.  Use  a  labeling  algorithm  to  determine  k(s,tlx.)  . 

w 

(k(s,tlx.)=l  if  s  and  t  are  connected 

J 

=0  otherwise. 

i.  Accumulate  results: 

S-S  +  k(s.tlXj)  . 

3.  Compute  final  estimate  for  K  replications: 
gK(s,t)  -  l-Fn_|C|(n,0)  +  H.S/K  . 


4.  Done. 


where  H  is  defined  in  step  2c  and 


atm  *  [Fm-rF|M|-l  *  7H7  <Fm‘Fm-l 

m  { 1  7i 

■  [Fm-rF|H|-i  +  <*->>  pi"o-p>n''”]/h 

bm  *  a«n  +  [Pmn-P)N'm]/H  .  (17 


Now  observe  that  one  can  write 


g(s.t)  =  l-FN|C|(N,p)  + 


r  rv^Vl  n N N— m 
l  l  *  P  (1-P) 

m=|M|  uJm 


Therefore, 

§K(s,t)  -  g(s,t)  *  1’fn_ | c |  (N.P)  +  £  .1  Ms.tlXj) 

J 

-  '  +  F„.|C|(N.P)  -  ”  f '  I..  Pm(1-p)N'm 


m= !M|  £CJ* 


s  H 


N-ICI  ,  K  »i 

1  I  l  iv  l  Ira  b  )  (O-P  0-P)  /HU  • 
m=  jM|  £eJ*  *  j*l  La£m’D£m;  J 
m 


Using  the  results  in  (7)  and  (9)  it  is  clear  that  this  bound  is  bounded 
above  by  2HQ/K  where  Q  is  a  function  of  the  number  of  nonconnected 


intervals  in  {[a  ,  b  ):  £cJ*  ;  m= |M| . . ,N- |C| }  .  In  particular, 

XIII  XIII  Ill 


the  worst  possible  arrangement  of  intervals  leads  to 


N-ICI 

Q  s  min{  V  R. 

n«|l 


,  T 


m-JMI  »’  ’  R"]>  ’ 


which  establishes  (1).  The  result  (ii)  follows  immediately  by  using 
(10)  instead  of  (7)  in  (15a). 


For  the  result  in  (iii)  note  that  steps  2a  and  2b  have  time 
complexities  0(n)  and  0(1)  ,  respectively,  per  replication.  Also, 
step  2c  can  be  performed  with  time  complexity  0(1)  per  replication 
using  the  cutpoint  sampling  method  described  in  Fishman  and  Moore  (1981). 
There  tables  need  to  be  computed  (with  time  complexity  0(n)  and 
space  complexity  0(2n))  prior  to  using  algorithm  C.l.  Step  2d  requires 
0(1)  time  per  replication.  Also,  the  identification 
of  operating  arcs  (steps  2e  and  2f)  takes  0(n)  time.  As  K-*«  ,  step 
2g  takes  0{ne)  time  on  average,  and  step  2i  takes  0(1)  time. 

Since  step  2h  based  on  a  depth-first  algorithm  takes  0(max(N,|V| ))  time 
the  result  in  (iii)  obtains.  This  completes  the  proof  of  Theorem  1. 

Observe  that  the  bounds  on  convergence  in  (i)  and  (ii)  are 
proportional  to  H»FN_ ^ C| (N,p)  -  F]M]1(N,p)  .  In  particular,  for 
fixed  N,  C  and  M  this  quantity  decreases  to  zero  as  p-*0  and  p-*l  . 

Also  if  C  and  M  are  unknown  one  can  substitute  0  for  iCI  ,  I  Ml 
or  both  in  algorithm  C.l,  which  continues  to  compute  correctly  but 
with  a  different  coefficient  in  the  bounds  in  (i)  and  (ii). 

The  quantity  Q  is  of  principal  interest  here.  It  depends  on 
the  arrangement  of  subintervals  In  the  unit  interval  assigned  to  the 
sample  outcomes  for  which  s  and  t  are  connected  for  each  m*IM|,...,N-ICI. 
If  for  given  m  these  subintervals  form  one  long  interval  with  no 
breaks,  then  Q*N-ICI-IMI+  1  .  However,  the  number  of  steps  necessary 
to  make  this  arrangement  prior  to  sampling  Is  comparable  to  the 
direct  computation  of  g(s,t)  by  total  enumeration.  Nevertheless, 
arrangements  that  are  clearly  preferable  to  others  hopefully  can  be  made  at 
relatively  small  cost.  These  assignments  are  a  topic  of  continuing 
research. 


2.  Stratified  Samp  line 


We  naturally  would  like  to  extend  the  techniques  of  Section  1  to 


the  estimation  of  g(s,t)  when  at  least  some  of  the  probabilities 


are  distinct.  This  extension  is  possible  if  one  exDloits 
certain  features  of  a  sampling  technique  known  as  stratified  sampling. 
Van  Slyke  and  Frank  (1972)  propose  this  technique  for  estimating  network 


reliability  characteristics  in  an  environment  in  which  q^,...,q^  are 


identical  and  samples  are  drawn  randomly  within  strata.  The  effect 


there  is  to  reduce  the  variances  of  the  resulting  estimators  below  the 


corresponding  variances  that  would  obtain  for  pure  random  sampling. 


This  section  retains  the  restriction  q^=...*q^=p  but  replaces  random 
sampling  within  strata  with  the  nonrandom  sequences  used  in  Section  1 


and  demonstrates  how  stratified  sampling  affects  the  bounds  on  error 


convergence.  Then  Section  3  shows  how  the  stratified  sampling  plan 


can  be  extended  to  the  case  of  distinct  q^,...,qN  . 


w  (s,t)  =  probability  that  s  and  t  are  connected  given 


that  exactly  m  arcs  are  operating. 


Now  observe  that  w  (s,t)  can  be  written  in  the  form 

m 


“»,<*•»>  *  -i,  I  k(s-tl2) 

'nr  m 


so  that 


g(s.t) 


l  %(*.*>  (Vo-p)N‘m 


N-ICI 


l-Fil.ir,(N,p)  +  l  pm(l-p)N"m  I  k(s,t|z) 


N-ICI 


m*|M| 


Algorithm  C.2  describes  a  procedure  for  computing 


l  K 

wjs.t)  =7  l  k(s,t|x.) 
m  K  j=1  J 


as  an  estimate  of  wm(s,t)  where  B  in  step  2c  gives  the  "number" 
of  the  combi  nation  of  A*m  operating  arcs  on  replication  j  . 
Algorithm  C.2 

Purpose:  To  compute  wm(s,t)  as  an  estimate  of  wm(s,t)  . 
Given:  N,  m,  Km,  s  and  t  (iMl^m^N-ICI). 

1.  Set  parameters:  n=N,  A=m  and  S=0  . 

2.  Perform  K  replications: 

m 

a.  For  1 n  set  x. .=0  . 

•  J 

b.  Select  v .  . 

J 

c.  Determine  the  "number"  of  the  combination  of  A  operating 
arcs  and  N-A  failed  arcs: 

b  •  i<;>  vaj  . 

d.  Determine  the  A-canonical  representation  a  ,  Bfc,...,B^  o 

using  (13).  If  B=0  ,  . 

e.  If  t>l  ,  then  for  v=l,...,*-l  ;  Bv=v-1  . 

f.  Set  operating  indicators: 

For  v=l,...,A  ;  set  i=B  +1  and  x.  .=1  . 

“  I J 

g.  Use  a  labeling  algorithm  to  determine  k(s,t|x.)  . 

\J 

h.  Accumulate  results: 

S*S+k(s,t|Xj)  . 

3.  Compute  final  estimate  for  K  replications: 

m 


ijs.t)  .  S/Km  . 


4.  Done. 


.“V.v.v  v  v  v  v ; 


N-|C| 


rNx_m,,  _xN-m 


gK(s,t)  =  i-fn_[C,(n,p)  +  m=IM  wm(s,t)(m)p  (l-p)‘ 


K=K IM |+ — +KN- |C | 


denote  an  estimate  of  g(s,t)  where  {w  (s,t);  m= |M| . N-|C|}  are 

computed  using  algorithm  C.2.  Then  Theorem  2  describes  bounds  on 
absolute  error. 

Theorem  2.  For  w  (s,t)  and  g,.(s,t)  algorithm  C.2  leads  to: 

"  m  i\ 


(i)  For  the  finite  sequence  {v.  =  -4^-  ;  j=l,...,K  } 

j  m 


'"m'5’1*  s2VKm 


(23a) 


where 


l5K(s,t)  -  9(s,t)I  s  2N'f '  ^  (!V<l-p)N‘l" 
K  m=|Mt  'Sn  m 


1  s  Qm  <  min[Rm(s,t)  ,  (JJ)  -  Rjs.t)]  m=|M|,...,N-|C| 


(23b) 


(11)  There  exist  infinite  sequences  {v.  ;  such  that 

J 

the  bounds  in  (23)  hold  with  c(log  Km) /Km  renlacinq  1/Km  . 
The  proof  of  (i)  relies  on  observing  that  for  m  operating  arcs 


l  m 

wJs.t)  s  v-  I  k(s,tlx.) 
01  m  j=l  1 


1 

l  4-  1  I r t-1  x  (V.) 

“j;  j-1  L(N,  •  (N,  >  i 

'nr  'm' 


and  on  using  methods  analogous  to  those  employed  in  the  proof  of 
Theorem  1  (i).  Part  (11)  follows  directly. 


Choosing  a  Samptz  Stze 


Note  that  for  a  given  set  of  sample  sires  {K  )  the  bound  in 

m 

(23b)  diminishes  to  zero  as  p-HD  and  p-*-l  .  However,  the  choice  of 
(K  )  clearly  affects  this  convergence.  Observe  that 
N-ICI 

°  '  nFLl5” 


where  Q  is  defined  in  (15b).  Therefore,  if  one  chooses 

Km  =  K(>m(l-p)N-m/H  m=IM| . N-lCl 

then  the  error  bounds  for  lgK(s,t)  -  g(s,t)l  in  Theorem  2  and 

A 

lgK(s,t)  -  g(s,t)l  in  Theorem  1  are  identical  for  the  finite  sequence. 
For  the  case  of  infinite  sequences,  the  dominant  term  in  the  bound  for 
lgK(s,t)  -  q(s,t) |  is  identical  to  the  bound  for  |gK(s,t)  -  g(s,t)|  . 
Moreover,  the  time  complexities  for  computing  gK(s,t)  and  gK(s,t) 
are  identical. 

An  alternative  approach  for  assigning  values  to  K  makes 

m 

considerably  better  use  of  a  priori  information,  when  it  is  available. 
Let 

X  *  mean  number  of  steps  needed  to  determine  whether  or  not 
s  and  t  are  connected  given  that  m  arcs  operate 
m- 1MI _ _ .N-ICI  . 

Then  for  a  specified  number  of  steps 


* 

K 


N-ICI 
j  K  A 

m=IMl  m  m 


the  assignment 


(24 


m*|Ml .... tN-|CI  . 


where 


«j  =  CQj(j)Pjn-p)N_j/^]1/2  J-1MI . H-ICI  , 

minimizes  the  bound  in  (23b).  Since  {Q  }  and  U  }  are  rarely 

mm 

known,  the  assignment  (25)  has  limited  practical  vlaue.  Note  that 
van  Slyke  and  Frank  chose  to  omit  consideration  of  the  computation 
time  when  choosing  their  assignment  of  { Km>  to  minimize  variance. 

There  is  an  additional  feature  of  stratified  sampling  that  makes 
algorithm  C.2  more  appealing  than  algorithm  C.l.  Observe  that  once 
(wm(s,t)  ;  m=lMl,...,N-lCl)  is  estimated  one  can  use  (22)  to  estimate 
g(s,t)  for  as  many  values  of  p  as  desired  without  the  need  for 
additional  sampling.  By  contrast,  algorithm  C.l  applies  for  only  one 
value  of  p  at  a  time.  Van  Slyke  and  Frank  also  mention  this  advantage 
for  stratified  sampling. 

3.  Unequal  Probabilities 

This  section  describes  how  stratified  sampling  can  be  extended  to 
the  estimation  of  g(s,t)  with  an  accelerated  convergence  rate  on  the 
error  bound  when  at  least  some  of  the  operating  probabilities  q^,...,qN 
are  distinct.  Diegert  and  Diegert  (1981)  also  describe  a  stratified 
sampling  plan  for  unequal  probabilities,  but  theirs  differs  significantly 
from  the  plan  proposed  here  both  in  terms  of  procedural  design  and 
sampling  mechanism.  In  particular,  their  plan  relies  on  independent 
replications  and  random  sampling  whereas  the  present  proposal  uses  the 
specialized  finite  and  Infinite  sampling  sequences,  discussed  in  Section  1, 
to  effect  the  accelerated  convergence  rates. 
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Let 

d  (s,t)  =  probability  that  s  and  t  are  connected  and 
m  that  exactly  m  arcs  operate. 

Observe  that 

z.  1-z. 

d  (s,t)  =  l  k(s,tlz)  n  Q,  O-qJ  (26) 

Z£Zm  ieP  1  1 

m 

so  that 

N  N-1 

g(s,t)  =  n  q.  +  l  d  (s,t)  .  (27) 

1=1  1  m= IMI  m 

Note  the  upper  limit  on  summation  N-1  replaces  N-ICI  to  account  for 

z.  1-z.  N 

the  probabilities  n  q - 1 ( 1 -q . )  for  N-|C|  <  I  z-  <  N-1  . 

i eP  1  i=l 

Algorithm  C.3  describes  a  procedure,  based  on  the  sampling 
plan  in  algorithm  C.2,  for  computing  dm(s,t)  as  an  estimate  of 
dm(s,t)  .  Then  our  estimate  of  g{s,t)  is 

gK(s,t)  =  n  qi  +  Tdm(s,t)  K=K(m|+...+Kn_1  .  (23) 

i  =  i  m= I  Ml 

N 

Recall  that  for  m  operating  arcs  there  are  (m)  possible 

N 

combinations  of  operating  arcs  and  that  Jm=(l . . ,(m) )  denotes  the 

set  of  all  indices  associated  with  these  combinations.  Also,  recall 
* 

that  denotes  the  subset  of  for  which  s  and  t  are  connected 

m  m 

when  m  arcs  operate.  Now  let  z^(t,m)  denote  the  status  (0  or  1) 
of  arc  1  on  combination  i  when  m  arcs  operate  i=l,...,N;  £eJm 
and  m=IMI ,. . . ,N-1  .  Then  Theorem  3  gives  the  relevant  bounds. 


Purpose:  To  compute  dm(s,t)  as  an  estimate  of 

dm(s,t)  *  probability  that  s  and  t  are  connected 
and  exactly  m  arcs  operate. 

Given:  N,  (q^ ;  i=l,...,N}  ,  m,  K  ,  s  and  t  .  (|M|  <  m  <  N-l). 

1.  Set  parameters:  n=N  ,  A=m  and  S=0. 

2.  Perform  K  replications: 

.  m 

a.  For  i=l,...,n  set  x.  .*0  . 

^  J 

b.  Select  v.  . 

J 

c.  Determine  the  "number"  of  the  combination  of  A  operating 
arcs  and  n-A  failed  arcs: 

B  *  L(J)  VjJ  • 

d.  Determine  the  A-canonical  representation  i  ,  B4,...,BA  of  B 

using  (13).  If  B=0  ,  £*A+1  . 

e.  If  t>l  ,  then  for  v=l,...,M;  Bv=v-1  . 

f.  Set  operating  indicators: 

For  v=l,...,A  ;  set  i=By+l  and  x^.*l  . 

g.  Use  a  labeling  algorithm  to  determine  k(s,tlx.)  . 

J 

h.  Accumulate  results: 

N  x. .  1 “X •  • 

S=S+k(s,t|x.)  n  q.  1 J (l-qi )  J  . 

3  i=l  1  1 

3.  Compute  final  estimate  for  K  replications: 


iswwwMW-Vi  W  TT rr  “'-.'V'--.' ;.-rvrl- l-  uVLP-  J7.r%,'-.,r  *'-.‘r-.^<L  •  %v"  •.*  oAw'.  <\*. 
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Theorem  3.  Consider  dm(s,t)  ,  based  on  algorithm  C.3,  and  gK(s,t) 
as  estimates  of  dm(s,t)  and  g(s,t)  respectively.  Then 

(i)  For  the  finite  sequence  (v.  =  ;  j=l,...,K  } 

J  m 


m 

2(U)  N  zAi,m) 

iam(S.t)  -  ys.t ^ 

111 


1-z. (£,m) 

(1-q,.)  (29) 


and 


N-J  (1!) 


N  z .  ( £,m)  l-z.(«.,m) 


Ig^s.t)  -  g(s,t)l  S  2  l  -g-  I  n  q1  (1-q. )  1 

K  m=  1 M I  *Sn  £eJ*  i=l  1  1 

m 

(ii)  There  exist  infinite  sequences  {v . ;  j=l,2,...}  for  which 

J 

(29)  and  (30)  hold  with  c(log  Km ) / Km  replacing  1/1^  . 
Proof.  For  (i)  note  that 


(30) 


N  z.(*,m)  1-z . (£,m) 

dm(s.t)  =  I  n  q,1  (1-q.) 

m  £€j*  i=l  1  1 

m 


(31) 


whereas 


,N,  Km 

.  m  -m 


N  x . . 


dm(s»t)  "  T~  l  k(s*tlxi)  n  9i1J  O-g* ) 
m  m  j=l  J  1=1  1  1 


u  „  1  m  N  z-(£,m) 

*  V  ^  t  l  JL)(yi)]  n  q«  (1-q,) 

lej*  J=1  L-fr-  ,  N  )  3  i=1  i  i 

'nr  'nr 


(32) 

l-z^.m) 


Therefore, 


N  z.(£,m)  l-z,(£,m) 

n  a  1  /i  ^  \  i 


iam(s.t)  -  <V,(s.t)|.  ig  Jj4  ^  q,’  (1-q,) 

m 

, 

'  r  l  tW  4-)<vi>  ■  ~tt  3' 
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(*% 

2V 

•Sn 


I 


AcJ* 

m 


N  z.(*,m)  l-z.(a,m) 

«  qf  '  (l-q,)  ' 

1*1  1  1 


* 


which  establishes  (29).  The  result  (30)  follows  by  summation  for 
m=  |M | , . . . ,N-1  .  Part  (ii)  follows  directly,  completing  the  proof. 

Note  that,  as  in  Theorem  2  and  3,  the  bound  (30)  diminishes  to 
zero  as  p*=max  q.  -*•  0  and  p*=min  q.  1  .  However,  in  contrast  to 

i  1  4  * 


the  case  of  equal  q^  ,  algorithm  C.3  requires  us  to  specify  the  unique 

<ll . ^  •  Therefore,  one  needs  to  use  algorithm  C.3  for  each  such 

set  of  probabilities  considered. 

There  remains  the  issue  of  how  to  choose  KjM| ,. . . ,KN_]  .  This  is 
not  a  trivial  problem,  as  Diegert  and  Diegert  show.  One  possibility  is 
to  choose 


k  (jj)(p«)"(i-p,)ll~m 


m=|M|,...,N-|C|  .  (32) 


Then  as  p*  -v  p*  the  bound  on  l§K(s,t)  -  g(s,t)|  in  Theorem  3 
converges  to  the  bound  on  J$K(s,t)  -  g(s,t)|  in  Theorem  1  for  the  finite 
sequence.  For  the  relevant  infinite  sequences  the  dominant  term  in  the 
bound  on  |gK(s,t)  -  g(s,t)|  converges  to  the  bound  on  |§K(s,t)  -  g(s,t)| 
proportional  to  (log  K)/K  .  Lastly,  the  time  complexity  for  computing 
9K(s,t)  converges  to  that  for  gK(s,t)  . 

While  (32)  enables  one  to  achieve  the  desired  accelerated  convergence 
as  in  the  case  of  equal  probabilities  in  Section  2,  better  assignments 
than  (32)  do  exist  and  conceivably  can  be  identified.  But  this 


identification  calls  for  additional  research. 

4.  An  Example 

Although  the  deterministic  upper  bounds  and  accelerated  convergence 
rates  make  apparent  the  superiority  of  the  proposed  method  over  any 
alternative  method  based  on  independent  replications  that  use  random 
sampling,  there  remain  insights  to  be  gained  from  an  illustration  of 
how  these  new  procedures  work  in  practice.  Accordingly,  this  section 
illustrates  the  performance  of  algorithm  C.l  using  the  van  der  Corput 
sequence  for  estimating  g(s,t)  in  the  network  of  Fig.  1  with  s=l  and 
t=20  .  Although  slower  in  convergence  than  the  sequence 

2i-l 

{Vj  *  ;  j=l , . . . ,K)  ,  the  van  der  Corput  sequence  enables  us  to 

proceed  sequentially  in  the  process  of  achieving  a  specified  accuracy. 


Insert  Fig.  1  about  here. 


The  network  in  Fig.  1  has  N=30  arcs  and  |V|=20  nodes.  The 
minimal  cardinality  among  all  minimal  cutsets  is  |C|~3  and  the 
minimal  cardinality  among  all  minimal  pathsets  is  |M|*5  .  Suppose  that 
all  arcs  have  operating  probability  p  and  that  the  objective  is  to 
estimate  g(s,t)  to  within  ±6  .  Then  for  the  van  der  Corput  sequence 

pr(|§K(s,t)  -  g(s,t)|  s{)=l  (33) 

for  all  K  2  K( 6 )  where 

K(fi)  -  min[j:  j«  -  2HQ(j  log2  j  +  1)  «  0;  j-2,3,...]  .  (34) 


Although  the  certitude  of  (33)  is  appealing  its  value  is,  at  the  moment. 


academic  since  Q  and  hence  K(6)  ,  are  unknown.  Therefore  we  use 
an  alternative  approach  for  studying  the  performance  of  algorithm  C.l. 

Consider  a  set  of  100  independent  experiments  on  the  Uh  of 

which  ( £«1 . 1 00)  the  van  der  Corput  sequence  (with  randomly  assigned 

starting  value)  is  used  with  algorithm  C.l  to  compute 
g|*^(s,t),...,g^(s,t)  as  estimators  of  g(s,t)  ,  where  K*  is 
a  large  integer.  Also,  suppose  that  on  experiment  *  the  quantity 

K^(«)  =  min(j :  |c|j^(s,t)  -  g(s,t)|  $6  j  <  i  <  K*)  (35) 

is  computed. 

Since  there  is  a  global  deterministic  upper  bound  K(6)  for  all 
van  der  Corput  sequences,  there  must  be  a  finite  number  of  replications 
associated  with  each  experiment  that  guarantees  the  accuracy  6  with 
certainty.  If  K*  is  sufficiently  large,  say,  K*  2  K( 6 )  ,  then 
KU)=KU)(fi)  i«l,...,100  are  the  required  numbers  of  replications  for 
these  100  experiments  and  a  study  of  the  properties  of  10  ,. . . 

should  provide  information  about  the  required  number  of  replications  to  achieve 
the  accuracy  6  with  an  arbitrary  van  der  Corput  subsequence  and  algorithm  C.l. 

For  our  experiments  we  used  p*.95  ,  leading  to 

g(l .20)  *  1-. 295414  *  10"  3  .  This  value  of  g(l,20)  is  actually  a 

22 

long  run  estimate  based  on  using  algorithm  C.l  with  K=2  *4194304 
points  generated  from  the  van  der  Corput  sequence.  Suppose  the  objective 

.  _c 

was  to  estimate  l-g(l,20)  to  within  -6  where  6*10  .  A  preliminary 

analysis  of  the  long  run  had  indicated  that  our  value  for  l-g(l,20)  was 
well  within  this  bound.  Moreover,  the  long  run  had  also  provided  evidence 


that  setting  K=2‘'1  =20971 52  in  (35)  was  reasonable. 

The  quantities  were  computed  on  100  independent 

experiments  each  using  the  van  der  Corput  sequence  with  a  randomly 

assigned  starting  value.  The  sample  mean  number  of  required  observations 

was  326469.58  with  estimated  standard  error  21289.28  .  The  smallest 

K ^  was  63343  ,  the  largest  was  857595  and  the  sample  coefficient 

of  variation  was  65.211.  Note  that  the  largest  value  was  well  within 
★ 

our  specified  K  .  Figure  2  shows  the  empirical  distribution  function 
of  the  data  with  the  deciles  marked. 


Insert  Fig.  2  about  here. 


As  a  basis  for  comparison,  suppose  we  had  used  algorithm  C.l  with 

pure  random  sampling  to  determine  v^,  v2 .  Then  Chebyschev’s 

inequality  gives  for  J  independent  replications 
prO^j  -  g|  $  6)  *  1  -  g(l-g)/J62 

where  §j=@j(s,t)  and  g=g(s,t)  .  If  5  is  sufficiently  small  so  that 
J  is  sufficiently  large,  then 


pr(lgj  -  9 Is  6)  s  2*(6  /J/g(l-gJ)  -  1 

where  *(•)  denotes  the  cumulative  distribution  function  of  the 
standardized  normal  distribution. 


Suppose  we  specify  a  confidence  level 

Pr( l9j  -  g I  s  6)  =  1  -  o  0  <  o  <  1  . 

Then  for  pure  random  sampling  the  approximating  normal  distribution 
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leads  to  a  minimal  required  sample  size 

J ( 6 ,ot)  =  Cd(a)/6]2  g(l-g)  (36) 

where 

d(a)  =  min[d:  4>(d)  =  1  -  a/2]  . 

For  a  confidence  level  1  -  a  =  .99  one  obtains 
J(. 00001,. 01)  s  19506035 

from  (36),  making  apparent  the  considerably  larger  sample  size  that 
random  sampling  requires  to  achieve  the  accuracy  ±1 0-^  99  out  of 
100  times  on  average. 

A  FORTRAN  program  based  on  algorithm  C.l  with  the  van  der  Corput 

sequence  used  an  average  of  1.391  microseconds  per  replication.  The  same 

program  using  a  sequence  of  independent  uniform  deviates  required  an 

average  of  1.329  microseconds  per  replication.  Therefore,  the 

van  der  Corput  sequence  requires  326469.58  x  .001391/60  =7.57  minutes 

to  achieve  the  desired  accuracy  whereas  independent  uniform  deviates 

require  19506035  x  .001329/60  *  432.06  minutes.  The  advantage  of 

using  the  van  der  Corput  sequence  for  estimating  g(s,t)  for  this  network 

is  apparent. 

5.  (s,T)  Connectedness 

We  now  turn  to  the  estimation  of 

g(s,T)  =  probability  that  node  s  is  connected  to  node 
t  for  all  •  teTCE-{s)  . 

On  replication  j  the  conditional  probability  of  (s,T)  connectedness 
Is 

Ms,T|x4)  «  n  k(s,t |x .)  .  (37) 
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Then  analogous  to  (5)  for  the  case  of  q^=...aq^=p  ,  this  suggests 
the  estimator 

3K(s,T)  =  1  I  k(s,T,xj)  ,  (38) 

3  1 

if  one  modifies  algorithm  C.l  as  follows: 

(i)  Replace  the  input  node  t  by  the  node  set  T  . 

(ii)  Replace  |C|  and  |M|  everywhere  by  unity. 

(iii )  Replace  k(s,t|x.)  by  k(s,T|x.)  in  steps  2h  and  2i . 

J  J 

The  motivation  for  modification  (ii)  is  to  account  for  the  fact  that 

there  may  be  different  sets  C(s,t)  and  M(s,t)  for  each  teT  and 

these  may  not  all  be  known.  If  they  are  known  one  can  replace  I  Cl 

by  min  |C(s,t)|  and  |M|  by  min  |M(s,t)|  . 
tcT  tcT 

Note  that  the  results  in  Theorem  1  continue  to  hold.  It  is  also 

possible  to  amend  algorithms  C.2  and  C.3  in  a  similar  manner  and  have 

the  corresponding  properties  in  Theorem  2  and  3  continue  to  hold. 
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